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Abstract 

It has previously been proven that finding the globally minimum energy configura- 
tion of an atomic cluster belongs in the class of NP-hard problems. However, this proof 
is limited only to homonuclear clusters. This paper presents a new proof which shows 
finding minimum energy configurations for heteronuclear clusters is also NP-hard. 



1 Introduction 

Atomic clusters are aggregates of atoms held together by the same forces that cause, for 
example, phase transition from vapor to liquid, formations of crystals, etc. Cluster sizes 
range from as few as three atoms up to several hundred atoms. The physical and chemical 
characteristics of a cluster often varies with its size. In fact, even the addition of a single 
atom can result in an entirely different structure. Only by successively adding more and 
more atoms will a crystal-like structure eventually be produced and some knowledge of the 
condensed phase attributes be determined 

The study of atomic clusters has steadily been increasing over the past decade 0. Of 
particular interest is the cluster conformation (structure) which has the lowest total internal 
energy E. Knowledge of this minimum energy conformation provides valuable clues relating 
to the chemical and physical properties of the cluster. Unfortunately, searching for the 
globally minimum energy state of a cluster has proven to be enormously difficult. Indeed, 
Wille and Vennik showed that locating the globally minimum energy state of a cluster 
of identical atoms — the homonuclear case — belongs in the class of NP-hard problems. This 
means there is little hope of exactly solving the problem in finite time for even moderate 
cluster sizes. 

The purpose of this paper is two-fold. First, it will be shown why existing homonuclear 
proofs, and work from other related problems, cannot be used for the heteronuclear case 
where not all of the atoms are identical. Secondly, a proof will be presented which does 
show solving the heteronuclear problem is NP-hard. 

2 Preliminaries 

The problem of finding this globally minimum energy structure is equivalent to optimizing 
E with respect to variations in all 3A^ — 6 degrees of freedom where is the cluster size. 
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One method of solving this optimization problem is to explore the potential energy surface 
(PES) composed of all possible cluster conformations. Each point on this surface represents 
a unique spatial arrangement of the constituent atoms. This multidimensional surface is 
characterized by numerous local minima, each indicating an energetically stable structure. 
If it were possible to enumerate all of these minima — and the saddles that link them — we 
could, in principle, describe the dynamics of chemical reactions governed by that surface. 
Unfortunately, enumerating all minima is extremely difficult because of their large number. 
Berry |^ indicates the number of geometrically distinct minima tends to grow exponentially 
with A^. Moreover, the number of permutational isomers grows factorially with A^. 

A number of general methods have been proposed for finding global minima on a PES, in 
particular, and on hypersurfaces in general. For example, Monte Carlo methods 0, eigen- 
vector following 0, evolution computation techniques lattice optimization/relaxation 
techniques |10|, and PES deformation techniques [11, 1^ have all been used with varying 



degrees of success. After formally defining the clustering problem, we will discuss some of 
these techniques in greater detail. 

The potential energy function for a cluster of A^ atoms is given by 

1 ^ 

y(r^) = 2 E - (1) 

hi 

where = (ri, r2, . . . , tn), Tj is the position vector of the ith atom and v{ri—rj) is a function 
representing the pairwise interaction between atoms. Lennard- Jones or Morse functions are 
frequently used for these functions. Our optimization problem of interest is the Discrete 
Cluster Problem (DCP) which is formally defined as follows:^] 

DISCRETE CLUSTER PROBLEM 

INSTANCE: Given finite number of points in real space, a distance d{i,j) G between 
two points an integer A^ G Z+ and a known potential energy function defined by (||). 
QUESTION: Is there a way to assign A^ atoms to A^ unique points so as to minimize the 
sum of their pairwise interactions? 

The potential energy of a cluster actually equals the sum of A^-body interactions. Re- 
stricting this sum to only two-body interactions provides only a qualitative approximation. 
Three-body interactions can be important, but are usually ignored in first order approxima- 
tions |jl3|. Although many-body interactions are needed for quantitative modeling, little is 
actually known about higher order terms. Many-body potential energy functions do exist 
though in practice, only two-body terms are used for the sake of computational speed. 



Consequently, the discussion here will likewise be restricted to the two-body case. 

It may appear that the large amount of work done with hard-sphere packing problems 
will be helpful in solving instances of DCP. Unfortunately, such is not the case because 
the objective of the two problems are quite different. Hard-sphere packings try to place A^ 
spheres in Euclidean space so that all can, without overlap, fit within as small a volume as 
possible [O. DCP deals with soft, compliant spheres which interact via pairwise interaction 



^The definition given covers both homonuclear and heteronuclear systems. 
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functions. Figure |1] shows a Lennard- Jones function, which is typical. Therein hes the major 
difference between hard-sphere and soft-sphere systems: the former has no preferred distance 
between the spheres while the latter does. Put another way, a hard-sphere packing algorithm 
attempts to minimize the interatomic distance r. Yet, a comparison with Figure clearly 
shows this does not yield the lowest energy state for an atomic pair. Hard-sphere packing 
studies can thus be expected to provide little help. 
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Figure 1: The scaled Lennard- Jones potential function v{r). The scaled interatomic dis- 
tance is denoted by r. 

An algorithm that searches for solutions to the homonuclear version of DCP was recently 
proposed by Hendrickson |1^. The general approach exploits structural information to 



decompose the global optimization problem into a set of smaller optimization problems. 
More specifically, objects in Euclidean space are represented as an undirected graph 
where the edge weights reflect the interobject distance. A subgraph is identified, the relative 
positions of the vertices are optimized, and the subgraph is then treated as a rigid body. 
These rigid bodies are ultimately combined to determine the overall structure. Such an 
approach won't work for the heteronuclear case, but this discussion will be differed until 
Section ^. 

Northby performed a lattice based search that takes an initial conformation and 
allows it to relax to an energy minimum. Essentially, a random search is conducted to 
compile a list of isomers. After pruning any geometrically equivalent isomers, the remaining 
conformations are relaxed under a Lennard- Jones pairwise interaction function. While this 
technique will work with heteronuclear clusters, it does depend upon a gradient search which 
makes it susceptible to stopping at local optima. 



Kostrowicki et al. ||TT| , presented a technique that "deforms" the PES in such a way that 
the number of local minima is dramatically reduced. This deformation requires solving a 
partial differential equation called a diffusion equation. All atoms are assumed to interact by 
a Lennard- Jones type of function which is approximated by a linear combination of Gaussian 
functions to permit expressing the answer analytically. Solving these diffusion equations is 
not necessarily easy as there are some pathological conditions. For example, numerical 
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problems can occur when the atoms are separated by large distances and so interact only 
weakly. This technique works best if one can suitably modify the boundary conditions for 
the diffusion equation, but care must be taken to ensure the minimum energy state is not 
removed from the deformed PES. 

Wales and Doye |T^ described another PES transformation technique called "basin hop- 
ping". Here the PES is converted into a set of basins of attraction for the local minima. 
This approach does has the advantage of not altering the energies of the minima. It has 
been used to find many of the local minima for Lennard- Jones cluster sizes up to = 110, 
The authors suggest some improvements that could be made to improve the efficiency of 
their approach. 

It is important to note that none of these techniques will, with certainty, guarantee a 
successful search for the globally minimum energy configuration. In the next section we will 
see why solving instances of DCP has proven to be so difficult. 



3 Complexity in Homonuclear Clusters 

An instance of DCP can be solved by exploring the PES associated the cluster, and any 
algorithm that does this can be used with both homonuclear and heteronuclear systems. A 
number of researchers indicate the number of local optima in a PES grows exponentially 



with A^ [g, |T^, |T8| . So, exactly how much effort is required to explore an exponentially large 
hypersurface? Some insight can be found from the theory of NP-completeness ||19|| . 

Suppose each time step a new point on the hypersurface is visited. Repeating this process 
until all points are visited — a procedure guaranteed to find the minimum energy state — will 
take an exponential number of time steps. Now consider an arbitrary NP-complete problem 
V. If we found an algorithm that could solve V in polynomial time, then this algorithm 
could solve every NP-complete problem in polynomial time. Unfortunately, no algorithm 
that solves an NP-complete problem has ever been found which executes in less than an 
exponential number of steps. This suggests the effort required to conduct a blind search 
of the hypersurface is similar to the effort required to solve an NP-complete problem. NP- 
complete problems are computationally intractable, which gives some idea of the difficulty 
one faces in solving DCP. (NOTE: This does not prove DCP is NP-complete. It merely 
shows the ramifications of having so many local optima.) 

NP-complete problems are defined as decision problems, that is, problems answered by 
either 'yes' or 'no'. NP-hard problems ask for the optimal solution to an NP-complete prob- 
lem. And, they have at least the same level of difficulty to solve as does the corresponding 
NP-complete problem. One particular NP-hard problem plays a pivotal role in the com- 
plexity analysis of DCP. It is the well known Traveling Salesman Problem (TSP) which is 
formally defined as follows: 

TRAVELING SALESMAN PROBLEM 

INSTANCE: A finite set C = {ci, C2, . . . , c^} of cities, and a distance d{ci, Cj) G for each 
pair of cities Cj, Cj G C. 
QUESTION: What permutation 

[C7r(l), C^(2), • • ■ , C,r(m)] 
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of C will minimize the tour length 

C,r(l)) ? 

Wille and Vennik used TSP to prove solving a homonuclear instance of DCP is NP- 
hard. It is worthwhile to examine their proof in some detail since our proof of NP-hardness 
for the heteronuclear case follows a similar line of reasoning. 

Instances of DCP can be expressed in graph-theoretical terms using a graph G = {V, E) 
with vertex set V and edge set E. The graph G is complete (i.e., the edge e = {u,v) G 
E W u,v G V). Furthermore, \V\ ^ A^. The edges are assigned weights w{e) W e E E 
where the weights reflect the interaction between vertices. An instance of DCP is therefore 
equivalent to selecting V G V with |\^'| = so that 




Figure 2: An example graph. The graph is dense and edges are weighted as given by Eq. 
The problem is to select A^ vertices so that the sum of the weighted edges is minimized. 

The proof of Wille and Vennik uses an undirected graph G = {V, E) with \ V\ = N{N—1) 
vertices (see Figure |^). Each vertex is labeled with (cj, Cj) where i, j = 1, . . . , N and i ^ j. 
This indicates city Cj is visited immediately after city Cj. Edges are then unordered pairs of 
the form 

e = {{ci,Cj){ck,ci)) 
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with edge weights 



w{e) 



-oo 



d{ck, ci] 
d{ci,Cj) 



if i^k,j 
if i = k,i 

if I = i, k = j 
if l = i^k^ j 
if k = jj^i 



(3) 



Selecting vertices, {(c7r(j), c,r(j+i)); i = ■ ■ ■ ,N, Ct,(^n+i) = C7r(i)}, so that the sum of the 
weights is minimal is equivalent to finding a minimal length tour (07^(1), 0^(2), • • • , c^(Af)) thus 
solving an instance of TSP. By restriction |jT9|, this also makes DCP NP-hard to solve. This 
completes the Wille and Vennik proof. 

The weight assignments given in Eq. (j^) require some clarification. In TSP, weight is 
equivalent to distance whereas in DCP weight is equivalent to pairwise potential energy. 
Any plausible solution to an instance of TSP can visit each city but one time. In valid tour 
moves, the edge weight equals the intercity distance. Disjointed tours have edges with zero 
weight — effectively removing that edge from the edge set E. This insures all tours can be 
defined as a permutation of m cities (see below). An illegal tour move has an edge with 
infinite weight. The total weight of the edges traversed in a tour measures the "goodness" 
of that tour; good tours have lower total edge weights. More specifically. 



1. i ^ k and j 7^ /. This defines a disjointed tour where the tour visits q followed by cj 
and Ck followed by q. However, it is not shown what other cities were visited between 
Cj and Ck- Therefore, it is not possible to describe the permutation and compute its 
tour length. 

2. i = k and j ^ l- This defines an illegal tour that leaves the same city to visit two 
different cities. 

3. j = I and i ^ k. This defines an illegal tour where two distinct cities visit the same 
next city. 

A. I = i and k = j. This defines an illegal cyclic tour among only two cities. 
5. l = i,k^j or k = j,l^i. These are legal tours. 



Each of the vertices in the minimal length tour is assigned a distinct atom from the 
cluster. In homonuclear clusters all one-to-one assignments of atoms to these N vertices 
are equivalent since the edge weights are based only on interatomic distance. This begins to 
explain why the Wille and Vennik proof and the Hendrickson algorithm |jl6| cannot be 
extended to the heteronuclear case. 

First, consider a homonuclear cluster with N > 2 atoms and suppose a new cluster is 
formed by swapping the spatial position of two atoms. It is not possible to tell any physical 
difference between the two clusters because all atoms are identical and the interatomic 
distances remain unchanged. Indeed, the new cluster will have a total energy identical to 
the original cluster because the pairwise interaction functions {v{ri — rj) in Eq. (|l|)) remain 
unchanged. Now suppose the cluster is composed of two distinct atom types, say Ar and Xe. 
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The repulsive and attractive forces experienced by an Ar-Ar atom pair differs from tliose 
experienced by Ar-Xe or Xe-Xe pairs even if all pairs are separated by the same interatomic 



distance |^ . Swapping atom positions in heteronuclear clusters changes the type of atoms 
which interact, altering the individual interaction functions, and giving a different total 
energy to the new cluster. 

With homonuclear clusters it is acceptable to consider atoms as simple identical spheres 
where only interatomic distances contribute to the total energy. It is natural to model this 
system as an undirected graph where the edge weights reflect forces derived solely from the 
interatomic distances. The search algorithm from [jl6l takes this approach thereby permitting 
the cluster to be optimized in stages by optimizing the relative positions in subgraphs. The 
complexity proof given in [Q] also made that assumption. In fact, the graph used for that 
proof was constructed specifically without requiring any pair type information to set the 
edge weights. That restriction was necessary to establish an equivalence between TSP and 
DCP. 

In heteronuclear systems, both atom type and distance determine pairwise forces so the 
corresponding graph must have edge weights that take both distance and atom type into 
consideration. Even the relative positions of vertices from a subgraph cannot be optimized 



without the weights being set in this manner. Consequently, search algorithms such as ||T6| 
cannot be used for heteronuclear clusters because the overall cluster structure depends on 
connecting rigid bodies, formed from optimized subgraphs, where the edge weights only 
consider interatomic distance. 

To apply the proof of Wille Vennik |^ to heteronuclear clusters, the graph would have to 
be augmented with additional vertices. For example, suppose in the original (homonuclear) 
graph two vertices i and j are connected by an edge e = Then w{e) is based solely 

on the interatomic distance between atoms i and j. The augmented graph would have 
a new vertex j' where the added edge e' = {i,j') has a weight w{e') computed from the 
interaction of two dissimilar atoms, spaced at a distance d{i,j') = d{i,j). However, now the 
mere selection of vertices — sufficient to solve a homonuclear DCP — does not guarantee 
the correct mixture of atom types present in the heteronuclear cluster. This restricts the 
proof from to only homonuclear systems. 



4 Complexity in Heteronuclear Clusters 

A different proof of complexity is needed for the heteronuclear clusters. This new proof 
makes use of the following NP-hard problem |jl9|: 

TRAVELING SALESMAN EXTENSION (TSE) 

INSTANCE: A finite set C = {ci, 02, ■ ■ ■ , Cm} of cities, a distance d{ci, Cj) G for each pair 
of cities Ci, Cj G C, a bound B G Z~^, and a "partial" tour 

= (C7r(l), C^(2), • • • , CTr(^K)) 

of K distinct cities from C, 1 < K < m. 
QUESTION: Can 9 be extended to a full tour 

(C7r(l), C,r(2)i • • • 5 C-^{K)i Cj^iK+l)! ■ ■ ■ , C^(m)) 
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having total length B or less? 

Consider a heteronuclear cluster which has a single atom of one type (denoted by a) 
and — 1 atoms of a different type (denoted by /?). A completely connected graph is 
G = {V,E) with \V\ = N{N — 1) is constructed with each vertex labeled as described in 
Section |[ Without loss in generality preassign an a-atom to all vertices with labels (ci, cj), 
j = 2, . . . , N and preassign a /5-atom to all remaining vertices. 

The search begins with a scan of all edges that touch vertices with a-atom assignments. 
Select the minimum weight edge.| This edge defines a partial tour (07^(1), c^(2))- Now select 
N — 2 more vertices, {(c7r(i), €7^(4+1)); i = 3, . . . , N, Cn{N+i) = Cn{i)}, so that the sum of these 
weights is minimal. (This search is limited to vertices with /?-atom assignments in order to 
maintain the proper mixture of atom types in the corresponding cluster.) Finding a minimal 
length full tour is thus equivalent to solving an instance of TSE. This proves that DCP is 
NP-hard to solve for heteronuclear clusters as well. 

5 Discussion 

The complexity in solving DCP forces researchers to use heuristic search algorithms. Hill- 
climbing algorithms are not expected to do well because the PES has many local optima 
and it is highly likely the algorithm will quickly stop at one of them. Conversely, stochastic 
search algorithms can be quite effective in such multimodal hypersurfaces. It is therefore 
natural to ask if any one particular algorithm stands out as doing especially well against 
DCP. 

This question is not easy to answer and we need to turn to the No Free Lunch (NFL) 
Theorem for an explanation. Essentially, this theorem says if some optimization (search) 
algorithm performs particularly well over a certain class of optimization problems, then it 
most likely will not perform as well over all remaining optimization problems. This means 
one cannot choose, for example, simulated annealing to use for DCP just because it happens 
to work well for, say, scheduling problems. Direct comparisons between algorithms are 
insightful. But, without a conscientious attempt to make the comparisons fair, the results 
may be inconclusive, or in the worst case, be completely wrong [^. Monte Carlo techniques 
have been dominant in the area of cluster studies, which is why newly proposed search 
algorithms are normally compared against them. 

Most notable among the new algorithms are the evolutionary algorithms, which conduct 
searches that mimic Darwinian evolution: a "population" of clusters evolve to a low energy 
state by altering existing configurations via stochastic reproduction operators. Natural se- 
lection determines which configurations survive to undergo further reproduction operations. 

Comparisons between evolutionary algorithms and Monte Carlo techniques have favored 
the former, although such comparisons sometimes lack sufficient mathematical rigor. For 
example, Zeiri 0| concluded that a genetic algorithm converges faster than simulated an- 
nealing after comparing the average of five runs — an unusually small sample size. Normally 
results should be averaged over a considerably higher number of runs to help remove any 
potential bias in the random number generators used in the algorithms. Nevertheless, there 

^As before, edges with weight are effectively removed from the graph. 
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is growing empirical evidence that says evolutionary algorithms consistently outperform the 
Monte Carlo techniques when applied against DCP |§||2^ I111"11H- 1^ appears as though 
evolutionary algorithms would be a good first choice search algorithm for cluster studies. 
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